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Second moment closure analysis 
of the backstep flow database 

By S. Parneix, D. Laurence 1 AND P. Durbin 2 


A second moment closure computation (SMC) is compared in detail with the di- 
rect numerical simulation (DNS) data of Le and Moin for the backstep flow at 
Re = 5,000 in an attempt to understand why the intensity of the backflow and, 
consequently, the friction coefficient in the recirculation bubble are severely under- 
estimated. The data show that this recirculation bubble is far from being laminar 
except in the very near wall layer. A novel ‘differential a priori 1 procedure was used, 
in which the full transport equation for one isolated component of the Reynolds 
stress tensor was solved using DNS data as input. Conclusions are then different 
from what would have been deduced by comparing a full simulation to a DNS. One 
cause of discrepancy was traced back to insufficient transfer of energy to the nor- 
mal stress by pressure strain, but was not cured. A significant finding, confirmed 
by the DNS data in the core region of a channel flow, is that the coefficient that 
controls destruction of dissipation, C €2 , should be decreased by a factor of 2 when 
production is vanishing. This is also the case in the recirculation bubble, and a new 
formulation has cured 25% of the backflow discrepancy. 


1. Introduction 

The flow over a backward-facing step has been probably the most popular sepa- 
rated flow test case of the past 20 years, for which numerous experiments (by Kim, 
Johnston, Eaton, Vogel, Durst, Driver, etc.) provide data on the effects of geome- 
try, inlet conditions, and Reynolds number. With the improvement of turbulence 
models and numerical methods, it is now generally possible to recover the reattach- 
ment length, but the intensity of the backflow and, as a consequence, the negative 
peak in skin friction are always underestimated by nearly a factor of 2 when second 
moment closures (SMC) are used. 

The recent DNS database of Le and Moin (1993) at Re = 5, 000, well corroborated 
by the experiments of Jovic & Driver (1995) and of Kasagi et al. (1995), is analyzed 
here to understand this severe defect common to all SMC. 

2. Full simulation of the backward- facing step 

The flow was computed using INS2D, a finite difference code in generalized co- 
ordinates written at NASA Ames Research Center. A fine, non-uniform grid of 
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120 x 120 cells was used to cover the region xjh — —3 to 35, x = 0 being the 
location of the sudden expansion and h being the step height. The inlet values for 
the mean velocities, Reynolds stresses, and dissipation were taken from the DNS 
database. The elliptic relaxation procedure of Durbin (1993) has been combined 
with the Speziale, Sarkar, Gatski (SSG) pressure strain model in the ‘neutral’ for- 
mulation as in Laurence et al (1995) (see Appendix). 



FIGURE 1. Streamlines (a) Second moment closure, (b) DNS 

Figure 1 shows the predicted streamlines compared to the DNS data. The reat- 
tachment is correct and a secondary bubble is found. However, the size of this corner 
bubble is much smaller in the simulation than in the data. In fact, if one looks at the 
predicted friction coefficient compared to the DNS and experimental data (Fig. 2), 
the intensity of the main recirculation is underpredicted by a factor 1/2 (the slight 
improvement shown by a dashed curved is discussed further in section 5). The 
stagnating flow between the two recirculations at x/h = 2 is also missed. Since it is 
believed that the underestimation of the secondary bubble is a consequence of the 
underestimation of the primary recirculation, we will concentrate in the following 
on curing the latter discrepancy. Note that the recovery after reattachment is also 
too slow — this is a problem in virtually all turbulence transport models. 




Figure 2. (a) Friction coefficient, o DNS (Le & Moin), ▲ experiment (Jovic 

Driver), SMC, modified SMC (cf. section 5), (b) DNS U-profiles at 

locations x/h — —3 (o ), 4 (□ ), 6 (o ), 10 (A), 15 ( <3 ) arid 19 (V), : model. 

The above observations are believed to reflect what can be expected from any 
state of the art SMC. In Fig. 2, the mean streamwise velocity U is shown. The center 
of the recirculation is well predicted, but its intensity is severely underestimated 
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(see station x/h = 4). Since Cf is too weak, this is not due to an overprediction 
of turbulent mixing in the near wall region, but rather to an underestimation of 
the entrainment from the shear layer. Since the flow in the upper layer splits at 
the stagnation point into the recirculation and the downstream flow, this same 
velocity defect is transported into the recovery region. Le & Moin noted, in good 
agreement with the Jovic & Driver experiment, that the log profile of the law of 
the wall was still not recovered at x/h = 19; this as a consequence of the low 
Reynolds number. In the Reynolds- averaged Navier-Stokes (RANS) computation, 
the recovery at x/h = 19 is, however, still slightly underestimated. 




FIGURE 3. (a) Zoom of U in recirculation at locations x/h = 1 (o ), 2 (□ ), 3 (o ), 

4 (A), 5 (<) and 6 (V), : model, (bed) Budget of U at locations (b) x/h = 2 

(c) x/h = 4 (d) x/h = 6, convection (• DNS, SMC) turbulent force (x DNS, 

SMC), viscous force (□ DNS, SMC), pressure force (* DNS, SMC). 

Figure 3 focuses on the recirculation and shows the budgets of the U momentum 
equation. The underestimation of the backflow is most severe at station x/h = 4. 
As expected, the adverse pressure gradient driving the backflow is fairly constant 
across the whole height and balances the viscous shear stress at the wall; hence, 
the pressure gradient determines the value of C/ min . Jovic and Driver (1995) found 
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that the minimum of Cf follows a iaminar like 1 law, C/ min = — 0.19i?e*~ i>/2 for 
Reynolds numbers between 5,000 and 50,000. It is very clear, however, that the 
recirculation is not laminar like, except for the very near wall region below the 
maximum of the reverse flow. 



FIGURE 4. (a) Zoom of V ill recirculation, legend cf. Fig. 3a. (bed) Budget of V 

at (b) x/h = 2 (c) x/h = 4 (d) x/h = 6, legend cf. Fig.3bcd. 

A possible explanation for the Reynolds number dependence, consistent with the 
fact that the turbulent Reynolds number remains high in the recirculation bubble, 
might be as follows. The pressure field is a consequence of the general form of the 
separated layer (which causes flow expansion and pressure rise). The form of the 
separated layer is determined by the turbulence. The Reynolds number influence 
noted by Jovic and Driver (1995) might come from a wider (compared to h) shear 
layer detaching from the step, causing a stronger adverse pressure gradient at lower 
values of Re . This view is supported by our observation that RANS computations 
were more sensitive to model changes in the shear layer than in the recirculating 
flow. The near-wall, viscous layer of reversed flow results from a balance of the 
pressure force and the viscous friction, and covers only a few percent of the bubble 
height— so this cannot lead to ‘laminar like 1 behavior. The turbulent Reynolds 
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k 2 e/v (i.e., ~ 10 v t /v) in the bubble is in the range 400 - 500 (Fig. 11), similar 
to its value in the shear layer. Incidently, this makes any low Reynolds ‘damping 
function’ ineffective in the main portion of the bubble. 

The momentum budget station at xjh = 2 shows that entrainment by the tur- 
bulent shear stress is underestimated right below the shear layer detaching from 
the corner. At xjh — 4 it was checked that uv and VdU/dy are by far the major 
contributors to the turbulent and convection terms. The shear is weakly opposing 
the recirculation (a pair of arrows in Figs. 3 and 4 indicates the position of U = 0) 
while advection is driving the recirculation in its upper part. Overall, the model 
seems only slightly to underestimate this turbulent shear force in the full simula- 
tion, even if at this stage the defect is traced to insufficient entrainment at the top 
part of the recirculation. However, we will see that uv itself is in error. 

The V component in the shear layer shown in Fig. 4 is severely underestimated. 
The turbulence force was found to be due almost entirely to v 2 . Again, right near 
the corner, at xjh = 2, the turbulent force is underestimated by about a factor of 2, 
although not far upstream at xjh = 0, the RANS results werejn accordance with the 
DNS data. At x/h = 4 and 6 (near reattachment), it is still v 2 which is driving the 
flow downward. Advection (inertia or streamline curvature) effects are negligible, 
showing that a reattaching flow is different from an impinging flow. Fig. 1 shows 
that the streamlines become smoothly tangent to the wall; some RANS simulations 
have produced a kink in the streamlines at this reattachment point (Hanjalic, 1996). 


3. Differential a priori tests 

S.l Reynolds stresses 

The Reynolds stresses will now be analyzed from two different sets of computa- 
tions. The first corresponds to the full simulation and explains the mean velocity 
budgets shown previously. The second is from a differential a priori test and per- 
mits an analysis of the true effects of the pressure-strain and transport models. The 
latter results are obtained by solving the full differential equations of each individ- 
ual stress UiU] one by one, while the other stresses and the mean flow are taken 
directly from the DNS database. 

An overall glance at Figs. 5-7 explains why a comparison using only the full 
simulation may entail erroneous conclusions: the full simulation could lead one to 
believe that the model has problems in the recirculation bubble, whereas the a 
priori test shows that discrepancies are mainly located in the shear layer. 

The u 2 streamwise fluctuation in the a priori test is overestimated in the shear 
layer, but improved in the recirculation when the correct mean velocity (used in 
the a priori test) enters its production. In both the full and a priori simulations, 
the v 2 component is seen to be underestimated in the shear layer, and elsewhere 
at station x/h = 4. On the other hand, the shear stress seems to be, on average, 
correct in the full simulation and overestimated in the a priori test. This is because 
in the former an erroneously small value of v 2 is entering its production term. The 
origin of the problem lies in insufficient return to isotropy in the SSG pressure-strain 
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FIGURE 5. t* 2 profiles (a) full computation, (b) a priori test, at locations x/h = 0.1 

(o ), 0.5 (□ ), 1 (o ), 2 (A), 4 (<) and 6 (V), : model. 




FIGURE 6. v 2 profiles (a) full computation, (b) a prion test, legend cf. Fig. 5. 


model, which should increase v 2 and decrease both u 2 and uv , but only in the shear 
layer. 

3.2 Budgets 

For the analysis of budgets, the DNS data has been processed in the same form 
as the elliptic relaxation model (see appendix); i.e., some anisotropy effects in the 
dissipation are lumped with the so-called pressure-strain term. 

For the budgets of u 2 in Fig. 8, the production terms coincide perfectly of course, 
since the mean velocities and Reynolds stresses other than u 2 , are taken from the 
DNS. That is the method of this differential, a priori test. The model for tur- 
bulent transport (Daly- Harlow) performs well, but the pressure gradient- velocity 
correlation (kf 12 ) is underestimated in the free shear layer. 

The budgets of v 2 (Fig. 9) show here again that in the shear layer the pressure 
correlation term is underestimated — at x/h = 2 by a factor of 2. The pressure term 
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FIGURE 7. uv profiles (a) full computation, (b) a priori test, legend cf. Fig. 5. 



FIGURE 8. A priori test: budget of u 2 at locations (a) x/h = 2, (b) x/h = 4, 
(c) x/h = 6, production (• DNS, SMC), [—eR t j/k] (x DNS, SMC), con- 
vection (■ DNS, SMC), transport (A triple correlations DNS, SMC), vis- 
cous diffusion (o DNS, SMC), pressure- deformat ion + dissipation anisotropy 

effects (V[n If 9 -ea+eRii/k} DNS, [kf tJ [ SMC) 
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FIGURE 9. 4 prion test: budget of e 2 at locations (a) x/h = 2, (b) x/ft = 4, 

(c) x jh — 6, legend cf. Fig. 8. 


includes pressure- transport effects (in the present case countergradient transport 
effects) which partially balance the turbulent diffusion terms, and which, as a con- 
sequence, are also underestimated by the model. Note that the production term is 
making a significant contribution, and since it is mainly composed of — v 2 DV/dy , 
underestimations of both v 2 and V (which is affected by v 2 ) self- amplify through 
this term. Near the wall, turbulent diffusion is generating the wall normal fluctua- 
tions, while the wall blocking effect is impeding them; the latter is represented by 
the elliptic relaxation effect (the homogeneous solution to Eq. 4 of the appendix is 
actually positive in this area). 

The budget of uv, on the other hand, shows an overestimation of the pressure- 
correlation, though again this compensates for an underestimation of the turbulent 
transport. At x/h — 4, near the wall, the production term is seen to change sign, but 
still Ttv remains approximately zero because of the strong transport term. Hence, 
in the narrow region between the maximum of the backflow and the wall, the mean 
flow is largely viscous, as seen previously. With increasing Reynolds number, one 
can expect the ratio of production to transport terms to become larger, and the 
turbulent shear stress would then decelerate this backflow, leading to a smaller 
peak in Cf. Because uv is countergradient with respect to the velocity gradient, 
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FIGURE 10. A priori test: budget of uv at locations (a) x/h = 2, (b) x/h — 4, 
(c) x/h = 6, legend cf. Fig. 8. 

the production of u 2 and k becomes negative in this area. 

4. Model parameters 

Fine tuning of models is often based on functions of the following parameters: 
anisotropy of the Reynolds stresses, A; turbulent Reynolds number, Ret : production 
over dissipation, P/e (sometimes the non-dimensional rate of strain, Sk/e is used); 
and turbulent lengthscale, k 3 ! 2 /e. In seeking improvements here, one should look 
for parameters that exhibit different values from those in simpler shear flows, for 
which the model should not be changed. The above parameters have been computed 
from DNS data to see if they are pertinent. 

The range of variation of A (Fig. 11) from 0.6 to 0.8 in the recirculation bubble 
shows no particularity; the Re t values in the range of 400 to 800 (Fig. 11) is too 
high to invoke low Reynolds effects; the ratio P/e (Fig. 12) decreases from 2 to 
1.5 in the shear layer, but is seen to be particularly weak in the lower half of 
the recirculation bubble. This last is a feature that is significantly different from 
near wall regions of boundary layers and should be considered further. Indeed, 
even negative production occurs along the wall from the reattachment to x/h = 4. 
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FIGURE 11. DNS profiles (a) Turbulent Reynolds number Re t = k 2 /(ve), at 
locations xjh = 2 (o ), 4 (□ ), 6 (o ), 8 (a ) and 10 (x ). (b) Anisotropy A — 
1 — 9 / 8 (A 2 — A 3 ), M and A 3 arc the second and third invariants of a j; , legend cf. 
Fig. 11a. 




FIGURE 12. DNS profiles (a) Production over dissipation P/e, legend cf. Fig. 11a. 
(b) Turbulent length scale L = 0.09 3 / 4 A ,3 / 2 /t , at locations x/h = 2 (o DNS, 

SMC), 4 (□ DNS, SMC), 6 (o DNS, SMC), 8 (a DNS, SMC) 

and 10 (x DNS, SMC), xy — . 


The production of dissipation is usually modeled as proportional to that of A*, and 
negative values might lead here to unphysical effects. 

5. Modeling dissipation 

Several attempts were made to increase the pressure- strain in the shear layer, 
but all resulted in a (sometimes dramatic) shortening of the reattachment length, 
without amplifying the strength of the recirculation. Though the previous analysis 
indicates that this is a route to pursue, the following only reports some success in 
improving the dissipation equation. 

An a priori test of the k — e equations was carried out by solving the coupled 
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Figure 13. A priori test: (a) k, (b) e, at locations x/h — 0.1 (o ), 0.5 (□ ), 1 (o ), 
2 (A), 4 (<) and 6 (V), : model, : modified model (e equation). 

k - e system with DNS data for uiUj and U. Dissipation is, of course, the exact 
source term for k, but k also has a strong relation to dissipation through the inverse 
timescale e/k in front of its source term: this is why the equations were solved as a 
coupled system. 

Figure 13 seems to indicate that k is overestimated and dissipation is correct in the 
coupled, differential test. In the recirculation bubble the effect of the source term in 
the e equation is destruction of e, since production is low. The ‘modification’ cited 
in the figure caption that is detailed below was intended to increase dissipation, but 
actually it leaves e unchanged and decreases k, bringing it closer to the DNS. This 
is because when the source term coefficient in e equation is decreased, the balance 
of the dissipation budget is re-established by a decrease of the time-scale, i.e. a 
decrease of k/e , that occurs by k decreasing with little change of t . 

Tuning of the dissipation equation has been a popular game for the past two 
decades, so it needs to be shown that the present modification should not deteriorate 
predictions in other flows, and what the rationale is behind it. 

Various procedures have been developed to enhance dissipation. In near wall 
flows below y + = 10, an extra viscous production term is usually included in low 
Re models. However, it is ineffective here because of the relatively high value 
of Re t . Another dissipation enhancement is the ‘Yap correction’ (Launder 1989), 
which consists in a positive source term in the dissipation equation that is activated 
whenever the turbulent lengthscale L is larger than the mixing length ny. Though 
rather ad hoc , this ‘Yap correction’ has been particularly effective for backstep 
or sudden expansion flows (Hanjalic 1996), and shows that something peculiar is 
happening to the dissipation that is still not understood. Indeed Fig. 12 shows that 
L is overestimated in the recirculation bubble near the wall. Although the Yap 
correction goes in the right direction, there is no justification for forcing L to be 
smaller than ny since the DNS data shows it to be considerably larger than this. 

Another way to increase the production of e in the near-wall region is to use 
the non-dimensionalized parameter P/e (Durbin 1993). This proved effective for 
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channel and boundary layer flow, but showed unfortunately high levels of numerical 
instability in more complex flows. Moreover, this does not suffice in the backflow 
region since P/e is fairly small. Durbin and Laurence (1996) recently proposed 
to replace this unstable term by a ratio of k and v 2 in the k — e — v 2 model: 

Ce x { 1 + a \ \J k/v 2 ) with = 1/30. In this study, we have generalized this idea 

in the full SMC by introducing the following: C f fl (l -F a\ y/Pk-e/\PsMc \ ) with 
ci\ = 0.035. Pk-e — 0.09 kTSijSji and Psmc = u x Uj S t 3 are respectively the k — e 
formula for production and the exact Reynolds-stress production. This correction 
has been found to have similar effects to P/e in the near- wall region of channel 
flow, without any numerical instabilities in more complex situations. It does not 
cure the underestimation of the backflow in the present case; of course, that was 
not its intent. 

Very little data is available concerning the dissipation equation budget aside 
from the channel flow DNS at CTR. From that data, the following adjustment to 
destruction of dissipation (which we called 'modified model 1 ) was devised: C f2 = 
1.83 */(/«) with I e = (P s + D £ )/2D e . This measures the weight of transport in the 
budget of e by using the imbalance of production minus destruction. The function / 
varies from 1 (for shear flows) to 0.5 when production is zero. In order to preserve 
numerical stability, we combined this modification with the same adjustment for 
C £} (multiplication by f(I e )). The following function / was chosen to avoid non- 
realistic coefficients: /(. r) = max(min(;r, 1 ), 0.5). The result for the backstep as 
concerns the C / profile was found to be modest (see Fig. 2), yet it is larger than 
any results obtained through modifications of the pressure strain model. 



FIGURE 14. Budget, of 5 for the Re r = 395 channel flow, all the terms have 

been multiplied by ( t/ + ) 2 , 'rapid part . 1 ( x [P\ + Pi + F 3 ], modelisation term: 

[1.44P*e/fc]) 'slow part 1 (+ [P.i — T], modelisation term: [— 1.83e 2 /t]), trans- 
port (□ DNS, modelisation term), viscous diffusion ), sum of the modeli- 
sation terms ( ). 
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The idea of re-adjusting C\ 7 actually came from the analysis of channel flow at 
Re r = 395, for which the budgets of dissipation are available (Mansour & Kim, 
private communication). The near wall region has been analyzed in detail by Rodi 
& Mansour, but what interests us here is the central part of the channel. It is well 
known to modelers that dissipation is underestimated in the core region, but with 
little consequence except that the modeled L is continuously increasing instead of 
leveling off just outside of the log-layer. Since € decreases as y“' 1 and the terms in 
its budget are as y“ 2 , all terms in this budget were multiplied by y 2 to produce 
Fig. 14. 

The dissipation budget, as discussed by Mansour, Kim & Moin 1988, or by Man- 
sour &; Moin 1993, comprises a viscous transport term, negligible in the core region, 
a turbulent transport, very well modeled here by gradient diffusion with a € = 1.3, 
and five source terms. It is known on fundamental grounds that these five terms 
cannot be clearly grouped into production and destruction terms. For the present, 
the terms involving gradients of velocity are grouped as Pi 4* P ‘2 + P* and com- 
pared to the ‘rapid 1 part of the model, e/k P k , while the remaining ‘slow 1 terms 
are compared to e 2 /k. The DNS values of k and e are used in the model terms, 
hence the jagged appearance of the model transport, due to double differentiation 
of this DNS data. It would seem from Fig. 14 that both constants C 5] and C f2 are 
severely overestimated; but again, the present split is arguable. It is, however, very 
clear that near the center of the channel, where the rapid terms go to zero, a value 
of C f2 — 1.83 is too large by a factor of 2. On the other hand, the transport is 
accurately modeled with the standard value cr £ = 1.3. 


6, Conclusion 

A detailed comparison of a SMC computation with the DNS data for the backstop 
flow at Re = 5,000 leads to the following conclusions: 

(1) The intensity of the backflow and the friction coefficient in the recirculation bub- 
ble are severely underestimated. 

(2) The recirculation bubble is far from being pseudo- laminar; an understanding of 
the problems encountered by SMC should, thus, be of general interest. 

(3) The SMC underestimates entrainment out of the recirculating bubble by the 
detaching shear layer. Thejnechanism is the following: pressure- strain <f > 2 2 gen- 
erates normal fluctuations v 2 , which create the transverse mean velocity V; this 
in turn provides a momentum impulse to the bubble. In the shear layer, < 5^221 e 2 , 
V are underestimated. 

(4) A new di fferential a priori procedure was used, in which the full transport equa- 
tions are solved one by one. 

(5) A modification was proposed to enhance dissipation in, and only in, the recircu- 
lation bubble. The new formulation cured 25% of the backflow discrepancy. 

(6) The model is supported by a significant finding from the DNS data in the core 
region of a channel flow: the constant related to destruction of dissipation, C 52 , 
should be decreased by a factor of 2 as production vanishes. 
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APPENDIX : Elliptic relaxation 

The Reynolds stress transport equation is written as: 

Dtuiiij = Pij + pij~ « 7 + Tij + vV 2 uJuj 


Pi] = ~UiU k dkUj - ujU k d k Ui 

Pij = -xt.djp - ujdip - (eij - uT Uj + ^ukdkpS,] 


Tij = -dk(u k uiu } + - UkpSij ) 


The term pij differs from the usual pressure-strain <j>ij since it includes a deviatoric 
dissipation tensor in the form 


Pij — 4*ij (^ £ ij 


The following neutral formulation for the elliptic relaxation is now obtained (Durbin 
and Laurence 1996): 

= ^ (4) 

k k k 

For homogeneous turbulence pij (= kfij ) in Eq. 4 reduces to pjj, for which any 
standard redistribution model <j>^j can be used. The SSG rapid model is 

= — C 2 dev (Pij) - C*dev{Dij) - C s kS tj (5) 


The coefficients are: 


^ 94 + 95 r, _ 94 “ 95 a 

C 2 = ; C7 3 - - 1 


C» = - g4 - 9s + ^ yfM 


The slow term is of the form 


= - [(Ci + 1 )an + C , ;dev(a lt « JL . J )] - 
Ci + 1 = - gi + fir*— ’ CJ = — — 


The dissipation equation is 


Cl,P-C t ,e 


Z?te = _£i_ 
The time scale, T, is defined as: 


+ 9t ^ + £^ )d , £ 


m /P ~v 
T = t/— + 36- 
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The length scale L appearing in Eq. 4 also is prevented from going to zero at the 
wall by using the Kolmogorov scale as a lower bound: 

[k? z / 3 / 2 

L = Ci YF + C ’?75- (9) 

Lastly, the Daly- Harlow expression for the turbulent diffusion is used: 

Tij = dt ( C', t ujxi m T d,„ nTUj ) (10) 

The constants used in this report are: 

C fl = 0.2, Of = 1.5, Cx = 0.1, C,, = 200, 

C et = 1.44, C*, = 1.83 

Also 

<i,j = dev(uiUj)/k, . 4 2 = a tj a,j, 

.4,-j = a,j «,*«*;,, .4 = 1- 9(.4 2 - A :i )/8 

was used in the text. 



